!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>


!     subroutine LIMBES_transport computes at the cell ni of the kinetic
!     equation v.grad(f) = vx df/dx + vy df/dy = vx(F_right - F_left)/dx +  vy(F_up - F_down)/dx 
!     The transport scheme is based on finite volume formulation.
!
!               F_up
!            ___________
!            |         |
!            |         |   
!     F_left |         | F_right  
!            |         |             
!            |_________|          
!              F_down
!     
!     These fluxes are computed at the cell ni by the subroutine F_flux(ni,F_right,F_left,F_up,F_down)
!     
!     LIMBES_vgradf(1:LIMBES_nmax, 1:LIMBES_nvmax)
!     -> v.grad(f) (real*8)
!     needed arguments
!     integer*8 ni cell index

subroutine LIMBES_transport(ni)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 ::  ni,vi,RLUD
  real*8 :: vx,vy
  real*8 :: F_right(1:LIMBES_nv)
  real*8 :: F_left(1:LIMBES_nv)
  real*8 :: F_up(1:LIMBES_nv)
  real*8 :: F_down(1:LIMBES_nv)


  !     computation fluxes on cell
  RLUD=1
  call F_flux(ni,F_right,RLUD)
  RLUD=2
  call F_flux(ni,F_left,RLUD)
  RLUD=3
  call F_flux(ni,F_up,RLUD)
  RLUD=4
  call F_flux(ni,F_down,RLUD)

  do vi=1,LIMBES_nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     LIMBES_vgradf(ni,vi)= vx*(F_right(vi) - F_left(vi))/SIMUL_dx + vy*(F_up(vi) - F_down(vi))/SIMUL_dx 
  enddo

end subroutine LIMBES_transport

!     --------------------------------------------------------------------
!     --------------------------------------------------------------------
!     --------------------------------------------------------------------

subroutine F_flux(ni,F_RLUD,RLUD)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: vi,nx,ny,ni,n1,n2,n3,nv,RLUD
  real*8 :: vx,vy,dt,dx
  real*8 :: psi,smooth,f1,f2,f3,CFL
  real*8 :: F_RLUD(1:LIMBES_nv)
  real*8 :: feq2,feq1

  nv=LIMBES_nv
  nx=SIMUL_nx
  ny=SIMUL_ny

  select case (LIMBES_transport_scheme)

  case("upwind1")
     !         LEFT/DOWN   RIGHT/UP
     !             |         |
     !             |    ni   |      
     !             |    n1   |
     !     for vx>0  F_right=f(n1) with n1=ni
     !     for vx<0  F_right=f(n1) with n1=LIMBES_neighbour_RLUD(ni,1)
     !     for vx>0  F_left=f(n1) with n1=LIMBES_neighbour_RLUD(ni,2)
     !     for vx<0  F_left=f(n1) with n1=ni
     !     extension straightforward for F_up and F_down
     do vi=1,LIMBES_nv
        vx=LIMBES_vx(vi)
        vy=LIMBES_vy(vi)
        !     F_right
        if (RLUD.EQ.1) then
           if (vx.GE.zero) n1=ni
           if (vx.LT.zero) n1=LIMBES_neighbour_RLUD(ni,1)
           f1=LIMBES_f(n1,vi)
           F_RLUD(vi) = f1
        endif
        !     --------------------------------------------------------------------
        !     F_left
        if (RLUD.EQ.2) then
           if (vx.GE.zero) n1=LIMBES_neighbour_RLUD(ni,2)
           if (vx.LT.zero) n1=ni
           f1=LIMBES_f(n1,vi)
           F_RLUD(vi) = f1
        endif
        !     --------------------------------------------------------------------
        !     F_up
        if (RLUD.EQ.3) then
           if (vy.GE.zero) n1=ni
           if (vy.LT.zero) n1=LIMBES_neighbour_RLUD(ni,3)
           f1=LIMBES_f(n1,vi)
           F_RLUD(vi) = f1
        endif
        !     --------------------------------------------------------------------
        !     F_down
        if (RLUD.EQ.4) then
           if (vy.GE.zero) n1=LIMBES_neighbour_RLUD(ni,4)
           if (vy.LT.zero) n1=ni
           f1=LIMBES_f(n1,vi)
           F_RLUD(vi) = f1
        endif
        !     --------------------------------------------------------------------          
     enddo
     !     --------------------------------------------------------------------
     !     --------------------------------------------------------------------
      case("upwind2")
!         LEFT/DOWN   RIGHT/UP
!             |         |
!             |    ni   |      
!       n2    |    n1   |
!     for vx>0  F_right=3/2*f(n1) - 1/2*f(n2) with n1=ni and n2=LIMBES_neighbour_RLUD(n1,2)
!     for vx<0  F_right=3/2*f(n1) - 1/2*f(n2) with n1=LIMBES_neighbour_RLUD(ni,1) and n2=LIMBES_neighbour_RLUD(n1,1)
!     for vx>0  F_left=3/2*f(n1) - 1/2*f(n2) with n1=LIMBES_neighbour_RLUD(ni,2) and n2=LIMBES_neighbour_RLUD(n1,2) 
!     for vx<0  F_left=3/2*f(n1) - 1/2*f(n2) with n1=ni and n2=LIMBES_neighbour_RLUD(n1,2)  
!     extension straightforward for F_up and F_down
         do vi=1,LIMBES_nv
            vx=LIMBES_vx(vi)
            vy=LIMBES_vy(vi)
!     F_right
            if (RLUD.EQ.1) then
               if (vx.GE.zero) n1=ni
               if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(n1,2)
               if (vx.LT.zero) n1=LIMBES_neighbour_RLUD(ni,1)
               if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(n1,1)
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               if (n2.EQ.ni) then ! if boundary, wall cell extrapolation
                  if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(ni,1)  
                  if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(ni,2)  
                  f2=LIMBES_f(n2,vi)
                  f2=two*f1 -f2
               endif
               F_RLUD(vi) = (three*f1 - f2)/two
            endif
!     --------------------------------------------------------------------
!     F_left
            if (RLUD.EQ.2) then
               if (vx.GE.zero) n1=LIMBES_neighbour_RLUD(ni,2)
               if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(n1,2)
               if (vx.LT.zero) n1=ni
               if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(n1,1)
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               if (n2.EQ.ni) then ! if boundary, wall cell extrapolation
                  if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(ni,1)  
                  if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(ni,2)  
                  f2=LIMBES_f(n2,vi)
                  f2=two*f1 -f2
               endif
               F_RLUD(vi) = (three*f1 - f2)/two
            endif
!     -------------------------------------------------------------------
!     F_up
            if (RLUD.EQ.3) then
               if (vy.GE.zero) n1=ni
               if (vy.GE.zero) n2=LIMBES_neighbour_RLUD(n1,4)
               if (vy.LT.zero) n1=LIMBES_neighbour_RLUD(ni,3)
               if (vy.LT.zero) n2=LIMBES_neighbour_RLUD(n1,3)
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               F_RLUD(vi) = (three*f1 - f2)/two
            endif
!     --------------------------------------------------------------------
!     F_down
            if (RLUD.EQ.4) then
               if (vy.GE.zero) n1=LIMBES_neighbour_RLUD(ni,4)
               if (vy.GE.zero) n2=LIMBES_neighbour_RLUD(n1,4)
               if (vy.LT.zero) n1=ni
               if (vy.LT.zero) n2=LIMBES_neighbour_RLUD(n1,3)    
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               F_RLUD(vi) = (three*f1 - f2)/two
            endif
!     --------------------------------------------------------------------
         enddo
     !     --------------------------------------------------------------------
     !     --------------------------------------------------------------------
  case("centered2")
     !         LEFT/DOWN   RIGHT/UP
     !             |         |
     !             |    ni   |      
     !             |    n1   |   n2
     !     for vx>0  F_right=1/2*f(n1) + 1/2*f(n2) with n1=ni and n2=LIMBES_neighbour_RLUD(n1,1)
     !     for vx<0  F_right=1/2*f(n1) + 1/2*f(n2) with n1=ni and n2=LIMBES_neighbour_RLUD(n1,1)
     !     for vx>0  F_left=1/2*f(n1) + 1/2*f(n2) with n1=ni and n2=LIMBES_neighbour_RLUD(n1,2) 
     !     for vx<0  F_left=1/2*f(n1) + 1/2*f(n2) with n1=ni and n2=LIMBES_neighbour_RLUD(n1,2)  
     !     extension straightforward for F_up and F_down
     do vi=1,LIMBES_nv
        vx=LIMBES_vx(vi)
        vy=LIMBES_vy(vi)
        !     F_right 
        if (RLUD.EQ.1) then
           n1=ni
           n2=LIMBES_neighbour_RLUD(n1,1)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,2)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 -f2
           endif
           F_RLUD(vi) =(f1 + f2)/two         
        endif
        !     --------------------------------------------------------------------
        !     F_left
        if (RLUD.EQ.2) then
           n1=ni
           n2=LIMBES_neighbour_RLUD(n1,2)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,1)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 - f2
           endif
           F_RLUD(vi) =(f1 + f2)/two
        endif
        !     --------------------------------------------------------------------
        !     F_up
        if (RLUD.EQ.3) then
           n1=ni
           n2=LIMBES_neighbour_RLUD(n1,3)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,4)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 -f2
           endif
           F_RLUD(vi) =(f1 + f2)/two
        endif
        !     --------------------------------------------------------------------
        !     F_down
        if (RLUD.EQ.4) then
           n1=ni
           n2=LIMBES_neighbour_RLUD(n1,4)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,3)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 - f2
           endif
           F_RLUD(vi) =(f1 + f2)/two
        endif
        !     --------------------------------------------------------------------
     enddo
     !     --------------------------------------------------------------------
     !     --------------------------------------------------------------------
  case("lax-wendroff")
     !         LEFT/DOWN   RIGHT/UP
     !             |         |
     !             |   ni    |      
     !             |   n1    |   n2
     !     for vx>0  F_right=f(n1) + 1/2*(1-CFL)*(f(n2)-f(n1)) with n1=ni and n2=LIMBES_neighbour_RLUD(ni,1)
     !     for vx<0  F_right=f(n1) + 1/2*(1-CFL)*(f(n2)-f(n1)) with n1=LIMBES_neighbour_RLUD(ni,1) and n2=ni
     !     for vx>0  F_left=f(n1) + 1/2*(1-CFL)*(f(n2)-f(n1)) with n1=LIMBES_neighbour_RLUD(ni,2) and n2=ni   
     !     for vx<0  F_left=f(n1) + 1/2*(1-CFL)*(f(n2)-f(n1)) with n1=ni and n2=LIMBES_neighbour_RLUD(ni,2)
     !     extension straightforward for F_up and F_down
     dt=SIMUL_dt
     dx=SIMUL_dx

     do vi=1,LIMBES_nv
        vx=LIMBES_vx(vi)
        vy=LIMBES_vy(vi)
        !     F_right  
        if (RLUD.EQ.1) then
           CFL=dabs(vx)*dt/dx
           if (vx.GE.zero) n1=ni
           if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(ni,1)            
           if (vx.LT.zero) n1=LIMBES_neighbour_RLUD(ni,1)
           if (vx.LT.zero) n2=ni
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vx.GE.0) then  ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,2)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 -f2
           endif
           F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)
        endif
        !     --------------------------------------------------------------------
        !     F_left
        if (RLUD.EQ.2) then
           CFL=dabs(vx)*dt/dx
           if (vx.GE.zero) n1=LIMBES_neighbour_RLUD(ni,2)
           if (vx.GE.zero) n2=ni
           if (vx.LT.zero) n1=ni
           if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(ni,2)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vx.LT.0) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,1)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 - f2
           endif
           F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)   
        endif

        !     --------------------------------------------------------------------
        !     F_up
        if (RLUD.EQ.3) then
           CFL=dabs(vy)*dt/dx
           if (vy.GE.zero) n1=ni
           if (vy.GE.zero) n2=LIMBES_neighbour_RLUD(ni,3)
           if (vy.LT.zero) n1=LIMBES_neighbour_RLUD(ni,3)
           if (vy.LT.zero) n2=ni
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vy.GE.0) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,4)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 -f2
           endif
           F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)
        endif
        !     --------------------------------------------------------------------
        !     F_down
        if (RLUD.EQ.4) then
           CFL=dabs(vy)*dt/dx
           if (vy.GE.zero) n1=LIMBES_neighbour_RLUD(ni,4)
           if (vy.GE.zero) n2=ni
           if (vy.LT.zero) n1=ni
           if (vy.LT.zero) n2=LIMBES_neighbour_RLUD(ni,4)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vy.LT.0) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,3)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 - f2
           endif
           F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)
        endif
        !     --------------------------------------------------------------------
     enddo
     !     --------------------------------------------------------------------
     !     --------------------------------------------------------------------
  case("lax-wendroff-P")
     !         LEFT/DOWN   RIGHT/UP
     !             |         |
     !             |   ni    |      
     !             |   n1    |   n2
     !     for vx>0  F_right=psi*f(n1) + 1/2*(1-psi)*(f(n2)+f(n1)) with n1=ni and n2=LIMBES_neighbour_RLUD(ni,1)

     !     extension straightforward for F_up and F_down
     dt=SIMUL_dt
     dx=SIMUL_dx

     do vi=1,LIMBES_nv
        vx=LIMBES_vx(vi)
        vy=LIMBES_vy(vi)
        !     F_right  
        if (RLUD.EQ.1) then 
           CFL=dabs(vx)*dt/dx
           if (vx.GE.zero) n1=ni
           if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(ni,1)            
           if (vx.LT.zero) n1=LIMBES_neighbour_RLUD(ni,1)
           if (vx.LT.zero) n2=ni
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vx.GE.0) then  ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,2)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 -f2
           endif
           feq1=LIMBES_feq(n1,vi)
           feq2=LIMBES_feq(n2,vi)
           smooth=(dabs((f1-feq1)/feq1) + dabs((f2-feq2)/feq2))/two
           if (smooth.LT.1.0d-2) then
              psi=CFL*smooth/1.0d-2
           else
              psi=CFL
           endif

           F_RLUD(vi) = psi*f1 + (one-psi)/two*(f2+f1)
        endif
        !     --------------------------------------------------------------------
        !     F_left
        if (RLUD.EQ.2) then
           CFL=dabs(vx)*dt/dx
           if (vx.GE.zero) n1=LIMBES_neighbour_RLUD(ni,2)
           if (vx.GE.zero) n2=ni
           if (vx.LT.zero) n1=ni
           if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(ni,2)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vx.LT.0) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,1)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 - f2
           endif
           feq1=LIMBES_feq(n1,vi)
           feq2=LIMBES_feq(n2,vi)
           smooth=(dabs((f1-feq1)/feq1) + dabs((f2-feq2)/feq2))/two
           if (smooth.LT.1.0d-2) then
              psi=CFL*smooth/1.0d-2
           else
              psi=CFL
           endif
           F_RLUD(vi) =  psi*f1 + (one-psi)/two*(f2+f1)
        endif

        !     --------------------------------------------------------------------
        !     F_up
        if (RLUD.EQ.3) then
           CFL=dabs(vy)*dt/dx
           if (vy.GE.zero) n1=ni
           if (vy.GE.zero) n2=LIMBES_neighbour_RLUD(ni,3)
           if (vy.LT.zero) n1=LIMBES_neighbour_RLUD(ni,3)
           if (vy.LT.zero) n2=ni
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vy.GE.0) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,4)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 -f2
           endif
           feq1=LIMBES_feq(n1,vi)
           feq2=LIMBES_feq(n2,vi)
           smooth=(dabs((f1-feq1)/feq1) + dabs((f2-feq2)/feq2))/two
           if (smooth.LT.1.0d-2) then
              psi=CFL*smooth/1.0d-2
           else
              psi=CFL
           endif
           F_RLUD(vi) = psi*f1 + (one-psi)/two*(f2+f1)
        endif
        !     --------------------------------------------------------------------
        !     F_down
        if (RLUD.EQ.4) then
           CFL=dabs(vy)*dt/dx
           if (vy.GE.zero) n1=LIMBES_neighbour_RLUD(ni,4)
           if (vy.GE.zero) n2=ni
           if (vy.LT.zero) n1=ni
           if (vy.LT.zero) n2=LIMBES_neighbour_RLUD(ni,4)
           f1=LIMBES_f(n1,vi)
           f2=LIMBES_f(n2,vi)
           if (n2.EQ.ni.AND.vy.LT.0) then ! if boundary, wall cell extrapolation
              n2=LIMBES_neighbour_RLUD(ni,3)    
              f2=LIMBES_f(n2,vi)
              f2=two*f1 - f2
           endif
           feq1=LIMBES_feq(n1,vi)
           feq2=LIMBES_feq(n2,vi)
           smooth=(dabs((f1-feq1)/feq1) + dabs((f2-feq2)/feq2))/two
           if (smooth.LT.1.0d-2) then
              psi=CFL*smooth/1.0d-2
           else
              psi=CFL
           endif
           F_RLUD(vi) =  psi*f1 + (one-psi)/two*(f2+f1)
        endif
        !     --------------------------------------------------------------------
     enddo
     !     --------------------------------------------------------------------
     !     --------------------------------------------------------------------
      case("mcd")
!         LEFT/DOWN   RIGHT/UP
!             |         |
!             |   ni    |      
!        n3   |   n1    |   n2
!     for vx>0  F_right=f(n1) - 1/2*(1-CFL)*(f(n2)-f(n1))*psi(smooth) with n1=ni and n2=LIMBES_neighbour_RLUD(ni,1)
!     for vx<0  F_right=f(n1) - 1/2*(1-CFL)*(f(n2)-f(n1))*psi(smooth) with n1=LIMBES_neighbour_RLUD(ni,1) and n2=ni
!     smooth=(f(n1)-f(n3))/(f(n2)-f(n1))   with :
!     n3=LIMBES_neighbour_RLUD(n1,2) for vx>0 and  n3=LIMBES_neighbour_RLUD(n1,1) for vx<0
!     for vx>0  F_left=f(n1) - 1/2*(1-CFL)*(f(n2)-f(n1))*psi(smooth) with n1=LIMBES_neighbour_RLUD(ni,2) and n2=ni   
!     for vx<0  F_left=f(n1) - 1/2*(1-CFL)*(f(n2)-f(n1))*psi(smooth) with n1=ni and n2=LIMBES_neighbour_RLUD(ni,2)
!     smooth=(f(n1)-f(n3))/(f(n2)-f(n1))   with:   
!     n3=LIMBES_neighbour_RLUD(n1,2) for vx>0 and  n3=LIMBES_neighbour_RLUD(n1,1) for vx<0
!     extension straightforward for F_up and F_down

!     psi=0  for smooth <= 0
!     psi=2*smooth  for 0<smooth<=1/3
!     psi=(1+smooth)/2  for 1/2<smooth<=3
!     psi=2  for  smooth>3
!     extension straightforward for F_up and F_down
         dt=SIMUL_dt
         dx=SIMUL_dx
         do vi=1,LIMBES_nv
            vx=LIMBES_vx(vi)
            vy=LIMBES_vy(vi)
!     F_right  
            if (RLUD.EQ.1) then
               CFL=dabs(vx)*dt/dx
               if (vx.GE.zero) n1=ni
               if (vx.GE.zero) n2=LIMBES_neighbour_RLUD(ni,1)
               if (vx.GE.zero) n3=LIMBES_neighbour_RLUD(n1,2)
               if (vx.LT.zero) n1=LIMBES_neighbour_RLUD(ni,1)
               if (vx.LT.zero) n2=ni
               if (vx.LT.zero) n3=LIMBES_neighbour_RLUD(n1,1)
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               f3=LIMBES_f(n3,vi)
               smooth=(f1-f3)/(f2-f1)
               if (smooth.LE.zero) psi=zero
               if (smooth.GT.zero.AND.smooth.LE.(one/three))  psi=two*smooth
               if (smooth.GT.(one/three).AND.smooth.LE.three) psi=(one+smooth)/two
               if (smooth.GT.three) psi=two
               F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)*psi
            endif
!     --------------------------------------------------------------------
!     F_left
            if (RLUD.EQ.2) then
               CFL=dabs(vx)*dt/dx
               if (vx.GE.zero) n1=LIMBES_neighbour_RLUD(ni,2)
               if (vx.GE.zero) n2=ni
               if (vx.GE.zero) n3=LIMBES_neighbour_RLUD(n1,2)
               if (vx.LT.zero) n1=ni
               if (vx.LT.zero) n2=LIMBES_neighbour_RLUD(ni,2)
               if (vx.LT.zero) n3=LIMBES_neighbour_RLUD(n1,1)                       
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               f3=LIMBES_f(n3,vi)
               smooth=(f1-f3)/(f2-f1)
               if (smooth.LE.zero) psi=zero
               if (smooth.GT.zero.AND.smooth.LE.(one/three)) psi=two*smooth
               if (smooth.GT.(one/three).AND.smooth.LE.three) psi=(one+smooth)/two
               if (smooth.GT.three) psi=two
               F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)*psi
            endif
!     --------------------------------------------------------------------
!     F_up
            if (RLUD.EQ.3) then
               CFL=dabs(vy)*dt/dx
               if (vy.GE.zero) n1=ni
               if (vy.GE.zero) n2=LIMBES_neighbour_RLUD(ni,3)
               if (vy.GE.zero) n3=LIMBES_neighbour_RLUD(n1,4)
               if (vy.LT.zero) n1=LIMBES_neighbour_RLUD(ni,3)
               if (vy.LT.zero) n2=ni
               if (vy.LT.zero) n3=LIMBES_neighbour_RLUD(n1,3)
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               f3=LIMBES_f(n3,vi)
               smooth=(f1-f3)/(f2-f1)
               if (smooth.LE.zero) psi=zero
               if (smooth.GT.zero.AND.smooth.LE.(one/three)) psi=two*smooth
               if (smooth.GT.(one/three).AND.smooth.LE.three) psi=(one+smooth)/two
               if (smooth.GT.three) psi=two
               F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)*psi     
            endif
!     --------------------------------------------------------------------
!     F_down
            if (RLUD.EQ.4) then
               CFL=dabs(vy)*dt/dx
               if (vy.GE.zero) n1=LIMBES_neighbour_RLUD(ni,4)
               if (vy.GE.zero) n2=ni
               if (vy.GE.zero) n3=LIMBES_neighbour_RLUD(n1,4)                        
               if (vy.LT.zero) n1=ni
               if (vy.LT.zero) n2=LIMBES_neighbour_RLUD(ni,4)
               if (vy.LT.zero) n3=LIMBES_neighbour_RLUD(n1,3)
               f1=LIMBES_f(n1,vi)
               f2=LIMBES_f(n2,vi)
               f3=LIMBES_f(n3,vi)
               smooth=(f1-f3)/(f2-f1)
               if (smooth.LE.zero) psi=zero
               if (smooth.GT.zero.AND.smooth.LE.(one/three)) psi=two*smooth
               if (smooth.GT.(one/three).AND.smooth.LE.three)  psi=(one+smooth)/two
               if (smooth.GT.three) psi=two
               F_RLUD(vi) = f1 + (one-CFL)/two*(f2-f1)*psi
            endif
!     --------------------------------------------------------------------
         enddo
     !     --------------------------------------------------------------------
     !     --------------------------------------------------------------------         
  case default
     print*,'error from LIMBES_transport: transport scheme not known'
     print*,'LIMBES_transport_scheme=',LIMBES_transport_scheme
     stop
  end select
end subroutine F_flux



